##########################################################################################

library(data.table)
library(optparse)
library(ArchR)
library(parallel)

##########################################################################################
option_list <- list(
    make_option(c("--comine_data_file"), type = "character"),
    make_option(c("--geneset_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){
    
    ## 整合atac和rna的文件
    comine_data_file <- "~/20231121_singleMuti/results/qc_atac_v3/germ_peak-gene/testis_combined_peak.combineRNA.qc.Rdata"

    geneset_file <- "~/20231121_singleMuti/results/celltype_plot/trajectory/positive/pct_0.25.list"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/celltype_plot/trajectory/positive/pct_0.25/smoothline"

}


###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

comine_data_file <- opt$comine_data_file
out_path <- opt$out_path
geneset_file <- opt$geneset_file

dir.create(out_path , recursive = T)
image_path <- out_path

###########################################################################################

a <- load(comine_data_file)
## testis_combined_peak_combineRNA
projHeme5 <- testis_combined_peak_combineRNA

dat_geneset <- data.frame(fread(geneset_file , header = F))

###########################################################################################
## 演化方向
cell_level <- c("SSC","Differenting&Differented SPG",
    "Leptotene","Zygotene","Patchytene","Diplotene",
    "Early stage of spermatids","Round&ElongateS.tids","Sperm"
    )

###########################################################################################
## 按照不同的时序分别去跑
## https://www.archrproject.com/bookdown/myeloid-trajectory-monocyte-differentiation.html
## 定义时序关系
trajectory <- cell_level

###########################################################################################
## 填补
projHeme5 <- addImputeWeights(projHeme5)

## 计算轨迹
## 发现每个单元格都有一个介于 0 到 100 之间的唯一伪时间值。我们排除具有NA值的单元格，因为这些单元格不是轨迹的一部分
projHeme5 <- addTrajectory(
    ArchRProj = projHeme5, 
    name = "Trajectory_time", 
    groupBy = "cell_type",
    trajectory = trajectory, 
    embedding = "UMAP", 
    force = TRUE
)

## 得到motif的名字
trajMM  <- getTrajectory(ArchRProj = projHeme5, name = "Trajectory_time", useMatrix = "MotifMatrix", log2Norm = FALSE)

###########################################################################################
## 画直线图
cpu <- 20

mclapply(dat_geneset$V1,function(gene){
  
    print(gene)
    motit_name <- grep( paste0( ":" , gene , "_") , names(trajMM) , value = T )[2]
    p1 <- plotTrajectory(projHeme5, trajectory = "Trajectory_time", colorBy = "GeneExpressionMatrix", name = gene, continuousSet = "horizonExtra")
    p2 <- plotTrajectory(projHeme5, trajectory = "Trajectory_time", colorBy = "MotifMatrix", name = motit_name, continuousSet = "solarExtra")

    image_name <- paste0( out_path , "/plotTrajectory." , gene , ".pdf")
    pdf(image_name)
    print(ggAlignPlots(p1[[2]], p2[[2]], type = "h"))
    dev.off()

},mc.cores=cpu)